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A method for computing inviscid hypersonic flow over complex configurations using unstructured meshes 
is presented. The unstructured grid solver uses an edge— based finite— volume formulation. Fluxes are 
computed using a flux vector splitting scheme that is capable of representing constant enthalpy solutions. 
Second— order accuracy in smooth flow regions is obtained by linearly reconstructing the solution, and 
stability near discontinuities is maintained by locally forcing the scheme to reduce to first-order accuracy. 
The implementation of the algorithm to parallel computers is described. Computations using the proposed 
method are presented for a sphere-cone configuration at Mach numbers of 5.25 and 10.6, and a complex 
hypersonic re-entry vehicle at Mach numbers of 4.5 and 9.8. Results are compared to experimental data 
and computations made with established structured grid methods. The use of the solver as a screening 
tool for rapid aerodynamic assessment of proposed vehicles is described. 


Introduction 

T HE current NASA/Industry effort toward next 
generation launch vehicles, namely, the Reusable 
Launch Vehicle (RLV) Program, has reinforced the 
need for rapid computational analysis of hypersonic 
re-entry vehicles. This program is being conducted 
under the NASA “better, cheaper, faster” philosophy. 
The initial design phase of the X-33 program required 
definition of a vehicle that could safely fly through all 
speed regimes and was less than 18 months in dura- 
tion. 

During the X-33 initial design phase, most of the 
aerodynamic and aeroheating data required for the 
design was obtained from wind tunnel tests, where 
models were rapidly tested across the subsonic to hy- 
personic speed regimes. The time to develop a suitable 
computational grid from an initial computer-aided de- 
sign (CAD) representation is typically several weeks, 
and small modifications to the geometry can require 
extensive re-working of the computational grids. As 
such, Euler and Navier-Stokes computational fluid dy- 
namics (CFD) capabilities are generally utilized at 
later stages of the design process; i.e., CFD codes are 
used to provide aerodynamics and heating information 
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on geometrically mature configurations. 

A rapid computational capability that provides a 
timely analysis of control surface loads and design 
modifications to the vehicle outer mold lines (OML) 
is needed to parallel wind tunnel efforts. Reason- 
able estimates of the aerodynamic characteristics for a 
flight vehicle are necessary in the early design stages, 
with the required accuracy increasing as the design 
matures. Unstructured mesh methods are ideally 
suited for rapid analysis as they have the capability 
to greatly reduce the time associated with grid gener- 
ation on complex configurations as compared to tra- 
ditional block structured grid methods. Unstructured 
mesh systems which provide automatic mesh gener- 
ation starting from a CAD definition and compute 
inviscid flow solutions have already been developed. 1-3 
These systems can handle complex configurations such 
as complete aircraft, and can produce reliable data 
using Euler flow models. However, most of the de- 
velopment work in the field of unstructured methods 
has been carried out for transonic or low supersonic 
flows. As these methods are extended to higher speed 
flows, they typically become less robust and require 
more dissipation and limiting to maintain stability. 

In this paper, the application of the unstructured 
mesh system FELISA 4 ’ 5 to the solution of flows about 
hypersonic re-entry vehicles is presented. The mathe- 
matical formulation of the solver is presented, and the 
issues unique to maintaining stable hypersonic compu- 
tations are discussed. The parallel implementation of 
the flow solution algorithm is also addressed. 

Three applications are presented to illustrate the 
developed capability. The first application compares 
FELISA results with inviscid LAURA 6 (a structured 
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grid code) calculations and experimental data for a 
sphere-cone configuration. In the second example, 
calculations are performed for a complex configura- 
tion and the resolution of flow features is compared 
with solutions calculated with viscous LAURA and 
inviscid DPLUR. 7 Lastly, the use of FELISA as a 
screening tool in preliminary design phases for rapid 
and accurate estimation of aerodynamic coefficients is 
demonstrated. Additional related applications have 
also been presented by Prabhu. 8 

Unstructured Mesh Generation 

The discretization of the three-dimensional compu- 
tational domain into tetrahedra is carried out by the 
mesh generator within the FELISA system. This re- 
quires that the geometry of the domain be defined in 
terms of an assembly of composite bi-cubic surface 
patches, which intersect along lines defined as com- 
posite cubic curves. For realistic geometries generated 
using CAD systems, the data files are sufficiently com- 
plex that tools for their efficient generation become 
necessary. 

The software package GridTool 9 allows the user to 
import geometries in a wide variety of formats, interac- 
tively perform the necessary geometry manipulations 
and output a FELISA compatible data file. These ge- 
ometry manipulations may include geometry clean-up, 
definition of intersection curves, re-parametrization of 
body surfaces, and definition of the far field domain 
boundaries. The grid density is controlled by means 
of a background mesh together with a system of point, 
line, and triangular sources. 

The first step in the mesh generation is the trian- 
gulation of the boundary surfaces using an advancing 
front method. The domain interior is then discretized 
using a constrained Delaunay algorithm which incor- 
porates exact integer arithmetic to resolve situations 
which are ambiguous when performed with standard 
floating point arithmetic. 10 

The use of the mesh generator in combination with 
GridTool has proven to be very reliable and easy to 
use. The speed at which elements are generated is 
above three million tetrahedra per hour using a sin- 
gle processor high-end workstation such as an SGI- 
R10000. The user-interactive part of mesh generation 
process is confined to defining the geometry and defin- 
ing the grid density. It is seen that the unstructured 
grid approach can give large reductions in grid gen- 
eration time as compared to structured grid methods, 
with typical times required for the generation of initial 
grids on the order of one to two days. The generation 
of multiple grids on the same or similar geometries can 
also be accomplished efficiently. 

Flow Solution 

The FELISA unstructured grid hypersonic flow 
solver uses a finite volume formulation. The ba- 


sic scheme is implemented using an efficient edge 
data structure 5, llj 12 and several flux vector/difference 
splitting options for computing the edge fluxes. The 
most successful and reliable option is the flux vector 
splitting proposed by Hanel et al., 13 which leads to a 
scheme that produces steady numerical solutions with 
constant total enthalpy (when the inflow conditions 
have constant total enthalpy). Second-order accuracy 
in smooth flow regions is obtained by linearly recon- 
structing the solution following MUSCL concepts. 14, 15 
Stability near discontinuities is maintained by locally 
forcing the scheme to reduce to first-order accuracy 
using the Local Extremum Diminishing criteria. 16 

Governing Equations 

The three-dimensional time-dependent Euler equa- 
tions written in integral form, over an arbitrary control 
volume D with boundary 3D, are expressed as 

4 f U dfi + f F k dS k = 0, (1) 

dt Jn Jd n 

where t denotes time, Su is the k - th Cartesian compo- 
nent of the outward area vector to 3D, and summation 
over repeated indices is assumed. The vector of un- 
knowns U and the vector of inviscid fluxes F p. are given 
by 
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( 2 ) 

The quantities p,p, «/., U, and H represent the density, 
pressure, k - th velocity component, total specific inter- 
nal energy, and enthalpy, respectively, and Sjk denotes 
the Kronecker delta. 

Spatial Discretization 

For each vertex i, a fixed control volume D a -, con- 
sisting of all tetrahedra sharing that vertex is defined. 
Considering equation (1) over each D a , and approxi- 
mating the integral over each triangular facet in 3D a 
by the value of the integrand at the centroid times the 
area we obtain 

if vdn+J2\(F l k + K + n)sl = o (3) 

J C2 i 7 G 

where denotes the set of triangular facets which 
form 3D a ; /, m and n, are the three nodes of a typical 
facet /; and = F^(U^) where denotes the value 
of the unknown vector U at node L By further ap- 
proximating the volume integral by U® times U®, the 
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volume of , we can write the following system of cou- 
pled ordinary differential equations for the evolution of 
the unknown vector at each node 

T T* 1 1 

~to + 777 £ 3 T + F ™ + n)Sl = 0. (4) 

In the above expression, the flux corresponding to 
a neighboring node appears as many times as facets 
share that node. For computational efficiency, it is 
convenient to re-arrange these contributions in such a 
way that instead of summing over the facets bound- 
ing fb', summation is performed over the mesh edges 
connecting i and its neighboring nodes. Thus, 

6 £ S i 

where Si denotes the set of mesh edges containing node 
«, e is a typical interior edge joining nodes i and l 
(Figure 1), and S e k is the k - th component of the area 
vector associated to edge e, which in turn is calculated 
as 

S e k = \j2 S l- ( 6 ) 

fer it 

In the above expression Tu is the subset of facets in 
that contain node L Noting that the sum of the 
area vectors over all edges belonging to Si must equal 
zero, expression (5) can be re-written for an interior 
node as 

r i\ T* 1 

^T + fS( F » + F')A* = 0 (7) 

6 £ S i 

where S e is the modulus of the area vector associated 
to edge e, = F*.S|, and S% — S%/S e . For interior 
edges, only one area vector needs to be stored for each 
edge. That is, the area vector required to form the 
equation for node l can be obtained from that required 
to form the equation for node i with a simple sign 
change. For boundary edges however, a special treat- 
ment is required to ensure that all the contributions 
defining the volumes associated with boundary nodes 
are considered. Details of this process are given in 
Peraire et ah, 17 and in practice it amounts to perform- 
ing an additional loop over all the boundary facets. 

Flux Vector Splitting 

The above scheme uses a centered approximation to 
the flux terms and is, therefore, non-dissipative. In 
order to produce a practical scheme which introduces 
sufficient damping to allow for strong shocks to be cap- 
tured, we follow a flux vector splitting approach and 
replace equation (7) by 
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Fig. 1 Set of triangular facets used to evaluate 
the area vector associated with an interior edge e 
joining vertices i and 1. 


f n ji i 

~tr + ^E (F t (ui) + F " (u ' ))5e = °- ( 8 ) 

6 £ S i 

Several choices for the numerical flux functions F^ 
are available in the literature. We have adopted the 
form proposed by Hanel et ah, 13 since our experience 
indicates that it provides one of the most competi- 
tive options when robustness, accuracy and cost are 
taken into account. In contrast, the unstructured 
mesh implementation of the flux difference splitting 
of Roe, 18 as used in the LAURA code, was found to 
be much less stable and in flows exhibiting strong nor- 
mal shocks, it often exhibited the so-called “carbuncle 
phenomenon.” 19 

The flux functions used here are 



±_f p(u n /c ± 1) 2 (2 =f U n /c), if |m„|<c 
^ \ p{u n ± \u n \)/2u n , otherwise 

( 11 ) 

respectively, where u n = and the sound speed 

c is given by c 2 = 7 p/p. An attractive feature of 
this form of flux vector splitting is that, whenever the 
enthalpy of incoming flow is uniform, the converged 
discrete solution also has constant enthalpy. 
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Fig. 2 Dummy vertices and elements used in the 
solution reconstruction process. 

Reconstruction 

The dissipation added using the above flux vector 
splitting is sufficient to make the overall scheme stable 
and allows for shocks of arbitrary strength to be cap- 
tured over just two or three cells. However, in smooth 
flow regions the scheme has poor accuracy, and re- 
quires unrealistically fine meshes to compute useful 
solutions. To obtain a scheme which retains good 
shock capturing properties, a geometric reconstruction 
scheme is employed. 

To illustrate this reconstruction procedure, consider 
(Figure 2 ) a typical interior edge e of length S e , joining 
nodes i and /, and define the dummy nodes i~ and 
by extending the edge e a distance S e at either side of 
nodes i and /, respectively. Let t{ denote the tetra- 
hedron containing node i, such that the solid angle 
defined by its three triangular facets containing node 
i includes the segment, or part of the segment, join- 
ing nodes i and i~ . In a similar manner, we define 
the tetrahedron ti such that it contains node l and in- 
cludes the segment, or at least part of the segment, 
joining nodes l and / + . The linear variations defined 
by the values of (/?, «/., H) at the four vertices of t{ and 
ti are used to determine the value of (p ) «/., H) at the 
dummy nodes i~ and , respectively. 

The modified algorithm, using reconstruction, is 
thus obtained by replacing ( 8 ) with 

r/TP 1 

+ ^E (F n (U+) + F n (U_)) ^ =°- (12) 

6 £ S i 

The unknown vectors U + and U _ are computed 
by first considering the vector of unknowns V = 
(/?, «i, 1 / 2 , U 3 , H) t and then evaluating U + = U(V + ) 
and U“ = U(V“), where V + and V are given as 
follows 

V+ = V 8 ' + |AV 
V- = V' - |AV 

AV = minmod(V l+ - V 1 , V 1 - V 8 ', V 8 ' - V 8 '”), 

(13) 


and where the minmod function returns the argument 
with smallest absolute value when all the arguments 
are of the same sign and zero otherwise. 

It is seen that by introducing the auxiliary vector 
V, the resulting limited scheme retains the ability to 
represent discrete constant enthalpy solutions. In ad- 
dition, it is observed that in smooth regions of the flow, 
the limited scheme can be interpreted as a central dif- 
ference approximation (7) with higher order artificial 
viscosity. 

Special care needs to be exercised near boundaries 
to appropriately choose the tetrahedron to be used in 
the determination of V at the dummy nodes. In this 
case a tetrahedron meeting the requirements specified 
above may not exist. A reasonable compromise that 
we have found to work well in practice is to choose 
the; 1 |fetrahedron that more closely meets these require- 
ments. 

Time Integration 

The equation set (12) is discretized in time using 
a forward Euler explicit time stepping scheme with 
a local time step determined from linearized stability 
analysis. For problems involving strong shocks, it is 
found that this form of determining the time step may 
lead to limit cycle oscillations that prevent the solu- 
tion from converging to machine precision. A remedy 
to that situation 20 consists of supplementing the lo- 
cal time step calculation with an additional check that 
enforces that the local monotonicity of the solution 
be preserved . 16 This additional check is found to be 
active for a very small number of points, but it is suf- 
ficeint to eliminate the limit cycle behavior. 

Boundary Conditions 

At the solid wall boundaries, the only contribution 
to the flux is through the pressure term. In addition, 
when the normal at the wall is well defined, the veloc- 
ity vector is forced to be tangent to the surface. Along 
trailing edges and other locations where a unique nor- 
mal is not defined, the velocity vector is not explicity 
modified. 

At the far field, boundary conditons are imposed 
using the flux splitting procedure described above and 
evaluating the appropriate flux function F± using far- 
held values. 

Parallel Implementation 

The above described explicit algorithm is easily par- 
allelized by decomposing the domain into subdomains 
which are then assigned to different processors and 
dealt with in a concurrent manner. The necessary 
communication between subdomains is accomplished 
using explicit message passing. The code can use ei- 
ther of the standard message passing libraries, PVM 21 
or MPI , 22 and runs on a number of platforms rang- 
ing from clusters of heterogeneous local workstations 
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to mainframe parallel computers such as the IBM SP2 
or the Cray T3D. 

Currently, the mesh generation is carried out in a se- 
rial mode. The mesh partitioning strategy is such that 
every edge in the mesh belongs to a single subdomain, 
whereas the nodal points on the subdomain boundaries 
are shared between two or more subdomains. The 
domain is decomposed using a simple coordinate bi- 
section algorithm. Within each subdomain, edges are 
divided into two groups. The first group is made up 
of those edges which have at least one node that re- 
quires communication (i.e., a node belonging to more 
than one subdomain). The second group is made up of 
the remaining edges. Typically, a loop over the edges 
in the first group is performed first. The values ac- 
cummulated at the boundary nodes are then broadcast 
to the neigboring subdomains. The interior edges are 
considered next, and then the information from the 
neigboring subdomains is received. This strategy al- 
lows for a considerable amount of communication and 
computation to take place concurrently, thus minimiz- 
ing latencies. The parallel efficiency of the code using 
16 nodes is higher than 95% for typical computations 
of the size reported in this paper. 

Results 

The current effort will demonstrate that the FE- 
LISA grid generation system and hypersonic Euler 
solver, FELISA_HYP, give reasonable and consistent 
predictions of aerodynamic performance and can be 
used effectively in rapid analysis of complex configu- 
rations. Predictions with FELISAJ4YP for a sphere- 
cone configuration and a preliminary RLV /X-33 lifting 
body configuration are compared with experimental 
data and with calculations made using a structured 
grid Navier-Stokes code, LAURA, 6 and an inviscid 
structured grid code, DPLUR. 7 LAURA has been ver- 
ified with Shuttle orbiter flight data 23,24 and a wide 
spectrum of configurations for supersonic to hyper- 
sonic flow conditions. 25 The parallel implementation 
of DPLUR has allowed it to be used as an aerodynamic 
screening tool for single-block configurations. 

The FELISA_HYP solver has both scalar and paral- 
lel implementations; however, all FELISA calculations 
presented herein were made using a single processor 
of an SGI-R10000 Onyx machine, as were the LAURA 
computations for the cone cases. The LAURA compu- 
tations on the RLV/X-33 lifting body were performed 
on the NAS Cray C-90. All DPLUR calculations were 
made on the IBM SP2 at LaRC. 

The strategy used in developing a block structured 
grid to perform a series of analyses on a given con- 
figuration using LAURA is to provide for the ’best’ 
grid that will be used. This means that each grid has 
spacing for viscous computations, and sufficient resolu- 
tion of the geometry to accurately capture the heating 
rates on the windside surfaces. LAURA utilizes an 


algorithm that adjusts the outer domain and spacing 
normal to the body to capture the shape of the outer 
bow shock for each case to be run, thus minimizing the 
number of complete grids that must be generated to 
analyze a particular configuration. For a typical grid 
on an X-33 configuration, the grid generation time is 4 
weeks (6 weeks were required on the first grid; lessons 
learned made subsequent discretizations faster.). 

The flexibility of the unstructured mesh approach 
offers a significant time advantage over the block struc- 
tured grid approach, and thus allows for a different 
strategy in developing meshes for series of computa- 
tions. The initial geometry definition for a FELISA 
mesh is done once for all grids, and produces a mesh 
that has the surface geometry adequately resolved, and 
an initial attempt at resolving the flow features of the 
first case. For typical X-33 geometries, this takes 
about 4 days; the first X-33 configuration that was 
analyzed required about 1.5 weeks. There are several 
strategies available for generation and adaption to a 
particular case. The ’brute force’ method was used in 
the examples in this paper; coarse mesh solutions are 
calculated, and the spacing distribution is revised with 
significant user intervention to capture the bow shock 
properly. Mesh enrichment, movement, and remeshing 
strategies have been developed for FELISA, but the 
results have not been satisfactory when the dominant 
flow feature is a strong bow shock. Solution adaptive 
strategies are under development which will allow for 
automatic resolution of the bow shock simultaneous 
with resolution of the flow features surrounding the 
body. 

Spherically Blunted Cone 

The first case presented is for flow over a 15° half- 
angle spherically blunted cone with a nose radius of 1.0 
inch. FELISA solutions were calculated at freestream 
Mach numbers of 5.25 and 10.6 and angles of attack 
of 0° and 15°. These results are compared to the ex- 
perimental results of Cleary 26 and computations made 
using DPLUR and LAURA in an inviscid mode. 

For the 0° cases, the FELISA solution was calcu- 
lated on a mesh of 267,459 points. The LAURA 
solution was calculated on a 61 x 33 axisymmetric 
structured grid. For the FELISA computations at 
Mqo = 5.25 and a = 15°, an initial mesh of 186,276 
points and a fine mesh of 693,910 points were used. 
For the M ^ — 10.6 and a = 15° case, the solution 
was calculated on a medium mesh of 328,527 points 
and the fine mesh used for the M ^ — 5.25, a = 15° 
case. The LAURA and DPLUR calculations utilized a 
61 x 19 x 33 grid, with the outer boundary adapted at 
each Mach number to the bow shock using the LAURA 
shock alignment routine. 

The theoretical value of C p at the stagnation point 
is 1.81 for Mqo = 5.25, and 1.83 for M^ = 10.6. All of 
the solutions presented here slightly underpredict the 
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Case 

Stag C p 
Moo = 5.25 

Stag C p 

Moo = 10.6 

Theory 

1.81 

1.83 

LAURA, 0° 

1.81 

1.82 

FELISA, 0° 

1.80 

1.82 

LAURA, 15° 

1.80 

1.82 

DPLUR, 15° 

1.79 

— 

FELISA, 15°initial 

1.78 

— 

FELISA, 15°medium 

— 

1.82 

FELISA, 15°fine 

1.80 

1.81 


Table 1 Stagnation C p values for 15° half— angle 
spherically blunted cone. 


stagnation point C p , with the LAURA results at 0° be- 
ing closest to the theoretical value. Table 1 shows the 
stagnation values of C p for all of the cases presented 
here. 

The comparisons of centerline pressure coefficients, 
shown in Figure 3, indicate good agreement for the 
computational cases. For the M ^ — 5.25 solutions, 
the largest differences in the C p values for the fine 
grid FELISA solutions as compared to the LAURA 
solutions are on the order of 2%, and are in the re- 
compression region of the windside of the cone at 15° 
angle of attack. In most regions, the difference is less 
than .5%. The solution on the initial mesh of 186K 
nodes gives a poorer comparison, due to inadequate 
resolution of the bow shock. Experience has shown 
that a mesh density of approximately 15 points is 
needed between the bowshock and the body to resolve 
the shock adequately; the initial FELISA mesh had 
approximately 10 points, while the fine mesh had ap- 
proximately 18. The comparisons for the M ^ — 10.6 
cases show differences between the LAURA and the 
FELISA solutions on the order of 6% in the recom- 
pression region, 2% elsewhere. The agreement (of all 
of the computations) with experiment is reasonable 
except for the 15° cases on the windside. 

RLV/X— 33 Lifting Body Configuration 

A preliminary Lockheed-Martin RLV/X-33 lifting 
body configuration is used in the remainder of the 
comparisons. The full vehicle is shown in Figure 4. 
It has twin vertical tails and fins; the engines are 
modeled by the box-shaped structure on the base. 
This configuration was evaluated during Phase I of 
the RLV/X-33 program with wind tunnel testing and 
computational analysis to predict aerodynamics and 
aeroheating for the full trajectory. The LAURA and 
DPLUR computations and the Mach 4.5 FELISA com- 
putations presented in this paper were run as part of 
the Aerothermodynamics Branch’s analysis program. 
The Mach 9.8 FELISA solution on this configuration 
was run specifically to compare with the LAURA and 
DPLUR results; the grid density was similar to a 
computation made with FELISA on an earlier con- 


O Experiment, Cleary 



a) Moo = 5.25 



b) Moo = 10.6 

Fig. 3 (^distributions for a 15° half— angle spheri- 
cally blunted cone. Experiment, LAURA, DPLUR, 
and FELISA 2HYP. 

figuration. 

RLV/X-33 Mach 9.8 Computations 

In this example, FELISA results are compared with 
viscous LAURA and inviscid DPLUR calculations for 
Mqo =9.8 and a = 40°. The base flow of the vehicle 
was not calculated. 

The computational grids for the FELISA and 
LAURA calculations are shown in Figure 5. The un- 
structured FELISA grid has 82K surface nodes, 716K 
volume nodes, and 4.2M tetrahedra. The LAURA vis- 
cous solution was computed on a 24-block structured 
grid system with a total of 1.85M grid points. The 
forebody portion of the grid was 65 x 65 x 65, and the 
12-block aftsection (including wings and vertical tail) 
was 101 x 248 x 65. The unusual shape of the exit 
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Fig. 4 Preliminary configuration of the Lockheed- 
Martin RLV/X-33. 

plane was generated using the shock alignment rou- 
tine within the LAURA code, and was copied for the 
FELISA.pase. The DPLUR data is for the forebody 
only, and used a single block, 51 x 65 x 41 grid. When 
compared to the structured computational grid, the 
FELISA mesh provided comparable resolution of the 
windside, much better resolution of the leeside of the 
vehicle, and somewhat coarser resolution on the wing. 
Also, the structured grid was a viscous grid and, as 
such, provided for finer normal spacing of grid points, 
and coarser axial spacing. 

Figure 6 shows the centerline pressure coefficients 
for the LAURA, DPLUR, and FELISA calculations. 
Overall, the agreement is very good, with small local 
differences in the recompression region of the forebody. 
Stagnation C p for all three cases is 1.81; thAlheoretical 
value is 1.83. 

A comparsion of the windside C p contours is given 
in Figure 7. The flow features of the FELISA solution 
compare very well to both the viscous LAURA and the 



a) FELISA mesh 



b) LAURA mesh 


Fig. 5 Computational grids for M oo = 9.8, 
a: z= 40° calculations on preliminary Lockheed- 
Martin RLV/X-33 configuration. 

inviscid DPLUR. The slight differences on the wing 
surface (Figure 7(a)) are attributable to grid density 
differences. The LAURA grid is more resolved in the 
wing region than this FELISA grid, due to the origi- 
nal purpose for each of the calculations. The LAURA 
calculation was made to resolve the heating levels on 
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Fig. 6 Cpdistribution on the symmetry plane of 
preliminary Lockheed— Martin RLV/X-33 configu- 
ration, Moo = 9.8 and (Y — 40°. 

the vehicle, while the purpose of the FELISA solutions 
was to generate first cut aerodynamic design data. 

The C p distribution on the leeside of the vehicle is 
shown in an expanded view (compared to the scale of 
Figure 6) in Figure 8. A significant variation between 
viscous LAURA and FELISA in the recompression re- 
gion is observed, due to the flow separation predicted 
by the viscous solution. The DPLUR inviscid and the 
FELISA compare very well up to an axial location 
of about 250. The FELISA grid is much denser in 
this region, and predicts a higher compression. The 
agreement with LAURA further down the body is 
completely fortuitous; the actual flow fields are very 
different due to the flow separation on the leeside. 

Comparisons of pressure coefficient contours on the 
leeside of the RLV/X-33 are shown in Figure 9. The 
difference in the crossflow shock structure between the 
FELISA inviscid and LAURA viscous solutions (Fig- 
ure 9(a)) is clearly evident. Comparison between the 
FELISA and DPLUR inviscid solutions (Figure 9(b)) 
shows similar crossflow shock structures. The DPLUR 
shock is slightly more smeared, and located further off 
the centerline of the body. This is due to the finer 
spacing of the FELISA mesh. In practice, structured 
grids are generated such that the leeside flow is not as 
resolved as the windside, so as to focus computational 
resources on the regions that have the most impact on 
the heating and aerodynamic forces on the body. 

Figure 10 shows pressure coefficient contours at a 
cut near the middle of the fuselage, comparing the 
FELISA inviscid results to the LAURA viscous and 
DPLUR inviscid. The LAURA solution clearly shows 
the separation as flow comes around the body. The in- 
viscid FELISA and DPLUR solutions remain attached, 
and generate a crossflow shock to turn the flow down- 
stream. The outer shock for the FELISA is more 
smeared; this is due to the the grid not being well 



a) FELISA and viscous LAURA solutions 



b) FELISA and DPLUR solutions 

Fig. 7 Cpcontours on the windside of prelim- 
inary Lockheed— Martin RLV/X-33 configuration, 
Moo = 9.8 and (Y — 40°. 
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Fig. 8 Leeside C p distribution on the symmetry 
plane of preliminary Lockheed— Martin RLV/X-33 
configuration, M oo = 9.8 and (Y — 40°. 

aligned with the shock and a coarse FELISA grid in 
the leeside bow shock region. Both of the structured 
grid solutions used LAURA’s shock alignment routine 
to capture the outer bow shock more efficiently, and 
it would be very difficult to generate an unstructured 
grid that provides as much alignment with the shock. 
Improvements to the adaptive capabilities of FELISA 
are needed to better align the bow shock, particularly 
on the windside of the vehicle. 

Figure 11 shows pressure coefficient contours for FE- 
LISA and LAURA viscous at a cut near the rear of 
the fuselage including the wing and vertical tail. The 
differences here are primarily in the region of the cross- 
flow shock and the outer bow shock. Qualitatively, the 
features compare very well, particularly in the flow 
around the wing. 

X -33 Aerodynamic Comparisons , Mach /.5 

The final example is a calculation of the aerody- 
namic forces over the complete X-33 configuration. 
The computations were performed at a Mach number 
of 4.5 in air, and at angles of attack of 0° and 10° . The 
experimental data shown for comparison is from a test 
performed in the LaRC Unitary Plan Wind Tunnel fa- 
cility, 27 and the data is courtesy of NASA LaRC. The 
test was run at M ^ — 4.50 and Re/ ft of 2 x 10 6 on a 
0.0105 scale model with transition artificially induced. 

The shaded surface and symmetry plane meshes 
used in the aerodynamic calculations are shown in Fig- 
ure 12. The FELISA grid used had a surface mesh of 
55K points, and a volume mesh containing 439K points 
and 2.56M tetrahedra. The same mesh was used for 
calculations at both angles of attack. The mesh ex- 
tended over the entire body, including the base region. 
The LAURA calculation used the same initial grid as 
the Mqo = 9-8 calculations and was adapted to form 
the grid shown in Figure 12(b) for the a = 10° case 
and a similar grid for the a = 0° case. Th|| base re- 



a) FELISA and viscous LAURA solutions 
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b) FELISA and DPLUR solutions 

Fig. 9 Cpcontours on the leeside of prelimi- 
nary Lockheed— Martin RLV/X-33 configuration, 
Moo = 9.8 and (Y — 40°. 
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a) FELISA and viscous LAURA solutions 


b) FELISA and DPLUR solutions 

Fig. 10 Cpcontours on forebody cut of prelim- 
inary Lockheed— Martin RLV/X-33 configuration, 
Moo = 9.8 and (Y — 40°. 


Fig. 11 Cpcontours on rear fuselage cut of the pre- 
liminary Lockheed— Martin RLV/X-33 configura- 
tion, Moo = 9.8 and cy = 40°, FELISA and LAURA 
viscous solutions. 

gion was not modeled in the LAURA computations. 

A plot of the normal force and moment coefficients 
vs. angle of attack is shown in Figure 13. The agree- 
ment of the FELISA normal force and moment co- 
efficients with the experimental data and the viscous 
LAURA calculations is very good. For the axial force 
coefficients, the agreement is slightly degraded. This 
is expected, since the FELISA calculations do not in- 
clude viscous effects. 

Conclusions 

A new algorithm for calculating inviscid hypersonic 
flow over re-entry vehicles using unstructured meshes 
has been developed. The flow solver uses a finite- 
volume formulation with an edge based data structure. 
Fluxes are calculated with the Hanel flux vector split- 
ting formulation, and a MUSCL gradient reconstruc- 
tion is implemented for higher-order accuracy. These 
algorithms for flux formulation provide more robust 
simulation capabilities for hypersonic flows as com- 
pared to earlier implementations in FELISA. The new 
algorithm has been implemented on both scalar and 
parallel architectures. 

The FELISA_HYP flow solver has been shown to 
be accurate in predicting both integrated aerodynamic 
characteristics and resolving flow features for complex 
re-entry vehicles. Examples of flow solutions over an 
RLV /X-33 configuration were presented. Comparisons 
of these solutions to experimental data and computa- 
tional solutions from established structured grid flow 
solvers provide the first, strong evidence for the valid- 
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ity of the FELISA_HYP algorithm in the code valida- 
tion process. The use of unstructured meshes permits 
faster mesh generation than traditional structured grid 
methods for complex configurations using state-of-the- 
art technology and expertise. 

The algorithm presented here enables rapid screen- 
ing for aerodynamics and pressure loads on space 
transportation vehicles in a time frame of two to four- 
teen days, including grid generation. The methodol- 
ogy presented in this work has been successfully ap- 
plied to preliminary phases of both the X-33 and X-34 
programs. The examples presented here are a small 
sampling of the work performed in the Aerothermo- 
dynamics Branch over the past two years. While the 
method described herein has been shown to be accu- 
rate and effective, there are several improvements in 
the grid generation and solution adaption areas that 
would further enhance the efficiency of the method. 



a) FELISA mesh, a = 0° and 10° 

Fig. 12 Body and symmetry plane grid for aero- 
dynamic calculations on preliminary Lockheed— 
Martin RLV/X-33 configuration, M oo = 4.5. 



b) LAURA grid, a = 10° case 
Fig. 12 Concluded. 



Fig. 13 Aerodynamic force coefficients for the 
Lockheed-Martin X-33 B1001A configuration at 
Moo = 4.5, (\ — 0° and 10°. 

These algorithm enhancements are the subject of on- 
going research. 
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